APPROXIMATION PROPERTIES OF THE g-SINE BASES 
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Abstract. For q > the eigenfunctions of the non-linear eigenvalue prob- 
lem associated to the one-dimensional g-Laplacian are known to form a Riesz 
basis of L 2 (0,1). We examine in this paper the approximation properties of 
this family of functions and its dual, and establish a non-orthogonal spec- 
tral method for the p-Poisson boundary value problem and its corresponding 
parabolic time evolution initial value problem with stochastic forcing. The 
principal objective of our analysis is the determination of optimal values of q 
for which the best approximation is achieved for a given p problem. 
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1. Introduction 

The p-Laplace operator or p-Laplacian, a generalization of the ordinary Laplace 
operator, arises naturally in applications from physics and engineering including: 
slow-fast diffusion related to particles [2] , superconductivity [5] , wavelet inpainting 
[23] , image processing [17] and game theory [22] ■ A typical application in the large 
p limit is a model for slow-fast diffusion for sandpiles [14] . Recently there has been a 
significant amount of research activity encompassing methods of approximation for 
solution of non- linear partial differential equations involving this operator [HI HI E] - 
The aim of the present paper is to further contribute to this activity by considering 
the particular case of the one-dimensional p-Laplacian and examine in detail the 
approximation properties of a generalized spectral method described as follows. 

Let p > 1. For z G K let [z] p_1 = z\z\ p ~ 2 . By extension from the linear 
case corresponding to p = 2, we define the one-dimensional p-Laplacian to be 
the differential operator A p u = Here u : [0,1] — > K is such that 

G i? 1 (0, 1). The corresponding p-Poisson boundary value problem is given 

by 

A v u(x) = g(x) < x < 1 

(1) P ~ ~ 
y ' u(Q) = u(l) = 

where g G £ 2 (0, 1). We also consider the related evolution equation 

r A , f „ , dW( X ,t) 

— — — = [A p u{x, t) - g(x)\ + v — 

( 2 ) u(x,0) = 0<x<l 
u(0,t) = u(l,t) = t>0 

that includes a stochastic forcing term where the noise intensity v > and W is a 
space-time Wiener process. Here 

(3) W(x,t) = J2l 3 n(t)Mx) 

nGN 

where /3 n are independent scalar Brownian motions, ip n are the eigenfunctions of 
the covariance operator Q and a n the corresponding eigenvalues, see [101 12T1 [T^] . 
For the case of space-time white noise we have Q — I. In the case of a deterministic 
system when v = 0, ([I]) is the steady state solution of 

Let q > 1. The so called q-sine functions are defined as the eigenfunctions of 
the g-Laplacian eigenvalue equation [131 023 HO] : 

-Mx) = ( g -1)A[ U ( 3 :)]'- 1 0<x<l 
u(0) = u(l) = 0. 

This family of functions and the corresponding problem Q was studied over 30 years 
ago by Elbert P~3] and later by Otani [5U], Bennewitz and Saito [7] in connexion 
with the computation of optimal constants in Sobolev-type inequalities. They gen- 
eralize in a natural fashion the 2-sine basis (corresponding to the linear case), and 
they have very similar periodicity and interlacing structures. 

In [8] analogues of the classical completeness and expansion theorems for the 
g-sine functions were established for q > 12/11. Specifically it was shown that they 
form a Riesz basis of L 2 (0, 1). This leads to the following question: what are the 
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approximation properties of this basis, as well as its dual basis, and how they relate 
to the approximation properties of the standard 2-sine basis? 

Below we address this question by examining approximation of the solutions of 
(JlJ and ([2| via projection methods with a g-sine and a dual g-sine basis, regarding 
q > 1 and p > 1 as free parameters. A main focus of attention is the determination 
of optimal values of q for which the highest order of convergence is achieved in 
a given p-problem. We demonstrate that standard properties of the 2-sine basis 
applied to the p — 2 problem (such as super-polynomial convergence when g(x) is 
smooth) are lost, when a g-sine basis for q ^ 2 is considered for ap^2 problem. 
As it turns out, the property of being a basis for the g-sine functions conceals a 
remarkably rich structure which is far from evident given the apparent simplicity 
of problem Q. 

Background material on the p-sme basis and its dual is considered in ££2j There 
we examine a matrix representation of the Schauder transform introduced in [5]. 
This will be crucial in our subsequent analysis as it gives rise to a stable procedure 
for constructing numerically both bases. 

In <|3] we find estimates for the approximation of square integrable functions in 
terms of their regularity. The dual g-sine basis turns out to have very similar ap- 
proximation properties as the q = 2 basis (lemma |2]). On the other hand, however, 
it is fairly simple to construct smooth functions such that their g-sine Fourier coef- 
ficients do not decay faster than a power —5/2 for q > 2 (lemma [3]). The latter is 
in stark contrast with the most elementary results in the numerical approximation 
of solutions of differential equations by orthogonal spectral methods. 

Section [4] is devoted to the p-Poisson boundary value problem. In theorem [5] we 
find explicit uniform bounds on the distance between any two solutions of ([I]), given 
the distance between the corresponding right hand sides. We then examine in detail 
the numerical computation of solutions of for source terms that are subject to 
various different regularity constraints. As it turns out, the estimates established in 
Theorem[5] appear to be sub-optimal. A more thorough investigation in this respect 
will be reported elsewhere. See [TT] for related results in the context of finite element 
approximation of the solutions of ([!]), including the higher dimensional case. 

In the final ^5] we study the numerical approximation of solutions to ([2| both in 
the deterministic and stochastic systems. We describe our discretization strategy 
and solve this problem for different values of p and v. Our results provide evidence 
on the performance of the g-sine basis for the solution of |2]), by showing the 
dependence on the parameter q of numerically computed L 2 residuals. 

2. The g-siNE basis and its dual 

The g-Laplacian eigenvalue problem Q, although non-linear, has a fairly simple 
structure. The eigenvalues are found to be A = {nir q ) q where ir q = qsi ^/ q ) ■ The 
first eigenfunction fi(x) associated to the first eigenvalue (ir q ) q is strictly increasing 
in [0, 1/2], decreasing in [1/2, 1] and it is even with respect to x = 1/2. It can be 
extended to an odd function (with respect to x = 0) in the interval [—1,1] and then 
to a 2-periodic C 1 function of R. If q > 2 then f"(x) is singular at x = 1/2. The 
eigenfunctions f n (x) associated to the eigenvalues (nir q ) q satisfy f n (x) = fi(nx) 
for all n > 2. 

In figure [I] we have plotted the first three eigenfunctions (top) and their corre- 
sponding derivatives (bottom) for (a) q = 1.4 and (b) q = 10. They typify the 
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case q < 2 (a) and q > 2 (b), respectively. For large q the basis functions /„ ap- 
proach zig-zag functions, which are the eigenfunctions of the oo-Laplace eigenvalue 
problem. 

Below we always assume that the family {f n }neN is normalized by the condition 
/4(0) = nn q and leave implicit the dependence of f n on q. In the special case q = 2 
we write e n (x) — \/2 sin(n7ra), so that {e„}„ e N is an orthonormal basis. 
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Figure 1. Approximation of the q-sine functions fj for j = 1,2,3 
(top) along with their derivatives (bottom) for (a) q = 1.4 and (b) 
q = 10. 



The Pythagorean identity generalizes to the q-sine functions [T3] as 
(5) [/i(ar)| fl + 7r-«]/i(a:)|« = 1. 

Integrating this differential expression for small enough x leads to the following 
explicit representation for the inverse function 

As we will see below, this representation plays a crucial role in the numerical esti- 
mation of f n (x). 

Let the Schauder transform, T q , be the linear extension of the mapping e„ i — > 
f n . Then T q : £ 2 (0, 1) — > L 2 (0, 1) is an invertible bounded operator for all q > 
12/11, [U Theorem 1]. Thus {f n }neN is a Riesz basis of L 2 (Q, 1) for such range of 
the parameter q. Further evidence presented in [5] suggests that in fact this is also 
the case for all q > 1, but at present this has not been proved rigorously. Unless 
otherwise specified we will assume from now on that q > 12/11. 

The property of a Riesz basis ensures that every g e L 2 (0, 1) is represented by 
a unique series expansion g = X^^Li a nfn which is convergent in norm. The q-sine 
Fourier coefficients, a n £ K, are given explicitly by a n = (g, /*) where {/^}neN is 



, 3 ,T*f*) for all 



the basis dual to {/„}„ g n- Since S jk = (fj,f%) = (T q ej,f%) = { 
j, k e N, then /* = (T~ 1 )*e n . It turns out that /* 7^ /„ for q ^ 2. 

The following matrix representation of T q is fundamental to our analysis. Let 



TqU) = AO') 



v2 / f\ (x) sin(j7ra;) da; 
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be the jth 2-sine Fourier coefficient of f\{x). Then the fcth 2-sine Fourier coefficient 
of f n (x) is given by 

fn(k) = V% / fi(nx) sin(fc7ra;) da; 
Jo 

(7) =2 fiifn) / sin(m7rna;) sin(fc7ra;) dx 

m odd ^° 

{T q (m) if mrt = fc for some m odd 
otherwise. 

Hence T q e n = Xim=i T ?( m ) e mn an d therefore T q has a lower triangular matrix 
representation in the orthonormal basis {e n } ne jq. See figure [2]- (a). 




Figure 2. In (a) we plot the distribution of the non-zero entries 
of a 1000 x 1000 truncation of T q . The insert corresponds to a 
20 x 20 truncation. The matrix entries are constant along each of 
the "quasi-diagonals" seen in the picture. In (b) we plot s against 
q where s is the numerically estimated H s regularity of the basis 
function fi{x). Note that at q = 2 the regularity is infinite. 



Remark 1. The basis of eigenvectors of the oo-Laplace eigenvalue problem are 
zig-zag functions, [HI Section 5]. In this case we can write ^(j) explicitly. As it 
turns out, 

g 

lim r q (j) = = Teoij). 

q— >oo 7 7T 

Let s > 0. Below we denote by iJp er [0, 1] the Sobolev space of 1-periodic func- 
tions g £ i 2 (0,l) such that £^(1 + j 2 ) s \g(j)\ 2 < oo. 



Lemma 1. Let f\ be the first q-sine function as defined above. If 1 < q < 2, then 
fl G H p 2 cr [0,l]. Ifq>2, then f t £ (J s <3/2 #per[0, 1] \ ff p 3 cr [0, 1] • Ifq> 4, then 
h $ H* el [0,l}. 

Proof. A straightforward argument involving integration by parts yields 



2y/~2 f 1 / 2 

Tqii) = — f"(x)sin(j<7rx)dx. 



3 k 
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Then 



(8) 



3 71 Jo 3 77 Jo 

_ 2y2 , 1/2 _ 2V2ir q 

„'2_2 Un x )\o — -2„2 • 



1/2 |/r(x)| 2 dx<l max \h(y)f = 4: 



Hence r q {j) = o(j~ 2 ) as j — > oo and /i 6 U s<3 ^ 2 iJp Cr [0, 1]. From |5| it follows 

that /('(z) - ft(/i(a;)) for h(y) = -v^V^l - y")^ . 
Let 1 < g < 2. Then 

1/2 

lf,"WI 2 rlT< ...... 

2 ye [0,1] 

so A ei^ er [o,i]. 

Let q > 2. Since lim^ii /"(x) = Too, then f x g i?p Or [0, 1]. Moreover, f x {x) > 
2x for < x < §. Then 

l/"W| 2 dx > j 1 ' 2 {2x) 2 ^ 1 \l - (2a;) 9 ) ^ dx. 

If g > 4, the integral on the right hand side diverges and so /i ^ -ffp CI .[0, 1]. □ 

Evidently the H s regularity for f\ (x) found in lemma [I] is not optimal. Fig- 
ure |2j-(b) shows a numerical estimation of the precise value of s(g), such that 
fi £ iJp Or [0,l] for r < s(q) and fi g Hp CI [0, 1] for r > s(q). The data for this 
graph was obtained by computing the decay rate of T q (j) for a large truncation of 
T q . A thorough investigation closely related to this lemma in the higher dimensional 
context can be found in [T^] and references therein. 

In the large q we have a limit of s = 3/2 and this is confirmed by Remark [I] At 
q = 2 we simply have f± = sin(7rx) and so f\ is in ifp Cr [0, 1] for all s. According 
to lemma [l] the curve should remain below s = 3 as q — > 2 + . For q > 3 the graph 
suggests /i ^ Hp C1 .[0,l]. However, if q < 3 it suggests s(q) > 2. As q — > 1 the 
regularity drops and the limit seems to approach s = 2. There is an interesting 
"peak" of regularity around q = 1.6 which we can not presently explain. 

The 2-sine Fourier coefficients of /^(x) are given by the n-th columns of (T^ 1 )*. 
This operator has an upper triangular representation in the 2-sine basis. Then 
/*(x) are trigonometric polynomials of order n. In fact, /|(x) = r ? (l) _1 sin(27rx) 
are parallel for all q > 1. Unlike the g-sine functions, not all dual g-sine functions 
have the same periodicity structure. 

We now describe a stable numerical procedure for computing these two bases. 
The g-sine functions can be approximated by first estimating /i using Numer- 
ical integration yields /-f 1 on [0, 1/2]. Although the integral is singular at y = 1, 
the value /^~ 1 (1) = 1/2 is known, so we do not need to consider quadrature points 
too close to this singularity. In our numerical procedure we chose a fine uniform 
grid and apply a cumulative Simpson's rule. This gives an approximation fx of 
/i for x £ [0,0.5] with a controlled tolerance and by symmetry we obtain f\ for 
x € [0, 1]. Note that /i is given on a non-uniform grid. The Pythagorean identity 
([5]) immediate yields the derivative /{ of f\ and hence we can use it to approximate 
the former with a /{, defined also on the non-uniform grid. 
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Once we have constructed /i and f[, we use periodicity and symmetry to find 
corresponding approximations of /„ and f n for n > 1 . This involves considering 
scaled copies /i and /{, to form /„ and j' n . Here the number of non-uniform grid 
points on [0, 1] grows with each n = 2 : N. 

To obtain /„ and f n on a uniform grid Xj — jh for j = : J, rather than the 
non-uniform grid that arises from the numerical integration, we have considered 
numerical interpolation by piecewise cubic polynomials. This gives {fn}n=i an( f 
{(fn) h }n=i defined at x — Xj. The former is the approximated basis and the latter 
the corresponding derivatives that we use for further computation. 

Figure [I] was generated with an implementation of the numerical scheme just 
described on an uniform grid with J = 4000 points (h — 2.5 x 10~ 4 ). The integral 
(|6| was approximated with 2 x 10 5 points. 

The dual basis is found from an N x N truncation, , of the Schauder trans- 
form. In practice, we first compute T q (j) and assemble . Then we define approx- 
imations (fn) h for n = 1 : iV, as the trigonometric polynomials whose kth 2-sine 
Fourier coefficients are the (n, k) entry of the matrix (T^) -1 . 

Remark 2. In our numerical approximation of the set of basis functions we are 
careful to fully resolve oscillations on the basis function fx, taking at least 20 mesh 
points per wavelength and for most computations 100 per wavelength. For N < 50 
we use h = (1007V)" 1 and so resolve each oscillation in the basis function with 100 
spatial points. For N > 50 we use h = (20N)~ 1 and so resolve with 20 points. 
We also examined convergence of orthogonality of the basis and dual in the spatial 
discretization and noted 0{h 2 ) for q « 10 through to 0(h 4 ) for q 1.4. 



3. Approximation of source terms 
Any given g £ L 2 (0, 1) can be approximated by either 

JV N 

9(x) « g* N (x) = ^2(g,fj)fj(x) or g(x) « g N (x) = ^(g, /*) fj(x) 



for large N. These two expansions converge as N — > oo in the norm of £ 2 (0, 1) and 
also pointwise for almost all x 6 [0, 1]. Unlike in the linear case corresponding to 
q = 2, the rates of decrease of || g — g^\\ and \\g — gj^W can be very different when 
q^2. 

Since the dual basis {/^} n eN comprises trigonometric polynomials, on the one 
hand we can formulate the following natural statement. 

Lemma 2. Let g e L 2 (Q, 1). For all JVeN, 



\9-9n\ 



< 



2V2 



N 



9(n)e n 



Proof. By definition g — g^ — X)^Ljv+i(^?*5' e «)(^ 1 g) le «- Since the matrix asso- 
ciated to T* is upper triangular (see Section [2l, 



(T*g,e n ) = J2[T*] nk g(k) = T q (k)g(nk) . 



k=l 
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According to Q, we have Y^k=i l T ?(^)l — ^73- Hence 

l(T;5,e„)| 2 = |5> g (fc)5(nfc)| = I £ t,(*0 1/2 t 9 (*) 1/2 ?(**)| 



k=l k=l 

oo 

< 



(E M fc )i) (E W*)iis(»*)i a ) < ^= E w*)ii£(»*)i s 

fc=l fe=l V fc=l 



Thus, 



OO OO 



n=JV+l v n =iV+l fe=l 



00 00 



K'w'vkY. E i^wii^i^ii^ii^Em^i E i^)i s 

2 fc=ln=JV+l Zy/Z k=l n=N+l 



2 00 

2 



<II^TA£k(*)l E I^HI^Tf E Wi\ 

V k=l n=JV+l 71=^+1 



□ 



Therefore the q-sine dual expansion of any g <= C°° converges super-polynomially 
fast. On the other hand, however, it is not difficult to construct examples of smooth 
functions g with a subsequence of q-sine Fourier coefficients decaying slowly. 

Lemma 3. Let g{x) = e±(x). If j is prime, then (g,f*) — T g (j)r g (l) . 

Proof. Since T q has an (infinite) lower triangular matrix representation in the or- 
thonormal basis {e n } n6 N, we can find the entries of the corresponding matrix rep- 
resentation of T^ 1 by pivoting and forward substitution (Gaussian elimination), ft 
is readily seen that T^ 1 is necessarily lower triangular and its diagonal should be 
constant and equal to r q (l)~ l . 

Assume that j is prime. According to ([7]), the only non-zero entries in the jth 
row of T q are r q (j) in the first position and r g (l) in the jth position. Therefore, 
the only non-zero entries in the jth row of (T g ) _1 are — in the first position 
and r 9 (f) _1 in the jth position. As the Fourier sine coefficients of f* are obtained 
from the jth column of (T*) _1 , the desired conclusion follows. □ 

By virtue of lemma [I] the prime g-sine Fourier coefficients of sin(7rx) for q > 2 
can not decrease faster than j~ 5 ^ 2 in the large j limit. Observe that this is in 
stark contrast with the most elementary results in the numerical approximation of 
solutions of differential equations by orthogonal spectral methods. 

Remark 3. The finite set of basis functions /„ and dual /* for n = 1 : N generate 
corresponding N dimensional subspaces Vat and of L 2 (0, 1). Instead of com- 
puting directly with these non-orthogonal bases, one can apply the Gram-Schmidt 
algorithm in order to obtain orthonormal bases of these subspaces. This has a 
numerical advantage of not needing to store both the basis and dual. We also 
considered this approach, however we found little advantage in terms of accuracy. 

Let us now consider various numerical tests on the approximation of regular 
functions by {/ n }neN and {/^} n eN- Once the bases have been obtained on uniform 
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mesh, we can examine their approximation properties numerically by looking at 
the decay of suitable residual for benchmark sources g € L 2 (0, 1). 

Figure [3] shows the typical outcomes of an experiment to determine the depen- 
dence in q of the L 2 (0, 1) residual. We have fixed here N = 40, varied q = 1 : 100, 
and computed residuals in the approximation of the following four functions: 

Sa(a0 = /i(a0 + (2.5)/io(a0 for q = 10, 
9b(x) 







x e [1/4,3/4] 
otherwise 



(9) 



=0*0 



(7/3)x x € [0,3/7] 

-{2l/4)x + 13/4 x e [3/7, 7/14] 

(21/4).t-2 x € [7/14,4/7] 

-a; + 11/7 x €[4/7,1] 

g d (x) = (q-l)(kn q rf 1 \f 1 \^ for q = 3, 

in (a)-(d) respectively. In the figure we include both the basis and its dual, as well 
analogous calculations with the orthogonalized bases of V40 and V^ . In (a) we see 
an optimal q op t = 10 as expected, in (b), (c) and (d) we increase the regularity of g 
and see g opt decrease with values of 4.25, 2.9, 2.55 for (b), (c) and (d) respectively. 

This experiment gives a general insight about the g-behavior of L 2 -residuals in 
the approximation of functions with different degrees of regularity by g-sine bases 
and their duals. For the simple functions considered, our calculations indicate the 
following general behavior of the g-sine basis: 

(i) as q — >• 1 the residual deteriorates, 

(ii) there is always a single minimum corresponding to an optimal q = (fopt, 

(iii) as q —¥ 00 the residual curve becomes asymptotically constant with no local 
maximum for q > 2. 

In contrast, the approximation error is almost constant in q for the g-sine dual basis. 
This is indeed a consequence of Lemma [2] and the fact that /* are trigonometric 
polynomials. Our tests indicate that there does not seem to be a clear advantage 
in orthogonalizing the basis or the dual basis for errors with an order of magnitude 
above 10 -5 . 



g{x) 

9b{x) 
9c{x) 
9d(x) 
sin(7rcc) 
Table 1. 



g = 1.8 
-0.4862 
-1.4266 
-1.9368 
-2.2778 



<7 = 2 (exact) 
-0.4905 (-0.5) 
-1.4474 (-1.5) 
-1.9952 (-2) 
NaN (-00) 



9 = 3 

-0.5005 

-1.4959 

-1.9826 

-1.9560 



9 = 5 

-0.5052 

-1.4842 

-1.5934 

-1.6237 



g = 10 

-0.5069 

-1.4438 

-1.4595 

-1.4988 



Rates of convergence in N for g with different degrees 



of regularity and selected values of q. 



In table[T]we have estimated a > such that ||<? — g^\\ < (3N~ a where /3 > is 
independent of N. The data indicates that for g £ Hp CI [0, 1], q pt ~ 2. Moreover, 
as q increases, we should expect a to decrease and stabilize always below 1.5 for 
g(x) = sin(7ra;). This is indeed suggested by lemma [I] figure [2]-(b) and lemma [3] 
and it is confirmed by the last row of the table. 

In figure[4]-(a) we have computed the residual \\g — <7io|| (red) and \\g — 520II (blue) 
for randomly generated g € ifp C1 .[0, 1] constrained to g(0) — g(l) = 0. To construct 
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(a) (b) 




Figure 3. We vary a g-sine basis and dual basis and examine how 
the residual changes using N — 40 modes for (a) g a combination 
of two 10-sine basis elements, (b) a piece-wise constant function 
<Zopt ~ 4.25, (c) a piece-wise linear continuous g pt = 2.9 and (d) 
a differentiable with discontinuous derivative function q opt — 2.55. 
See ^. 



a random function g such that the i/p OI .[0, 1] norm is finite, we find g~(j) — ajfij 
where a,j = (1 + j 2 )~ s ^ 2 \j\~^~ s for some small 5 > and f3j are independent 
identically distributed N(0, 1). Two realizations of functions g obtained in this way 
are shown in the inserts in figure [4j Note that q opt is achieved at different places 
but close to 2. In figure |l]-(b) we examine q op t over 200 realizations of functions in 
ifp O1 .[0, 1]. For any fixed realization and value of N, q opt ^ 2 although in the limit 
this optimal parameter appears to be close to 2. For N = 100 we also show the 
results of 1000 realizations. In this case, the mean value is closer to 2 although the 
variation is still large. In figure |4]-(c) we show the mean values where functions are 
taken in H poT [0, 1] for s = —0.5, 0, 0.5, 1 and 2 with 200 realizations. For s < 2 the 
variability in q opt is far larger. For fixed N, q op t > 2. 

The observed outcome of this experiment strongly support the conjecture that 
<fo P t typically approaches 2 as the regularity of g is increased. 
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(a) (b) (c) 




Figure 4. In (a) 1? residual as a function of q where a g-sine 
basis is used to approximate two different randomly generated g £ 
Hp CI [0, 1] satisfying g(0) — g(l) = 0. Note that the optimal g op t 
occurs at different values. In (b) and (c) we examine q opt as N 
increases. We include the mean and standard deviation over 200 
realizations of random functions subject to the same constraint. 
For N — 100 we also show the results from 1000 realizations. 



4. Numerical solution of the p-Poisson equation 

We now address the question of approximating the solutions of ([I]) by means 
of a g-sine and dual basis. In view of lemmas [2] and [3j we begin by determining 
uniform estimates on how sensitive this solution is under perturbations of the right 
hand side. Analogous questions have certainly been considered in more general 
frameworks, however here we focus on the explicit calculation of the constants 
involved. 

A key ingredient in the estimates presented below is the fact that can be 
integrated explicitly. Let the Volterra operator 

Vg(x) = [ g(t) dt. 
Jo 

Note that Vg £ H 1 ^, 1) for all g £ L x (0, 1). Furthermore V : L s {0, 1) — > L*(0, 1) 
is a contraction operator for all 1 < s, t < oo and its norm can be explicitly 
determined [71 Theorem 1.1]. Let 

Jo 

Then hgfa) is a continuous function, decreasing in 7, for all fixed g £ L (0, 1). Let 
min Vg(x) < 70(g) < max Vg(x) 

x£[0,l] x£[0,l] 

be the unique root such that /i s (7o(s0) = 0. Then 

(10) u(x) = J lVg(r) - 70(5)]^ dr = V (lVg{r) - 7 o(<?)]^) (*) 

is the unique solution of ([!]). 

We firstly establish concrete Holder estimates on h g {^). Without further mention 
in this section we will fix r = and denote by || • || s the norm of L s (0, 1). In the 
case s = 2 we will continue suppressing the sub-index. 
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Lemma 4. Let g <E L 1 (Q, 1) and — \\g\\i < 7 < [i < ||.9||i- TTien 

(M-7)< l f^[h g (j) - h g (^)} 0<r<l 

(M-7) r <2'- 1 [^(7)-^(/i)] r>l. 

Proof. Suppose first that < r < 1. From the graph of Jz] r for \z\ < M it is 
readily seen that [z] r - (wf > rM r - x (z - w) for all -M <w < z < M. Then 

(/i - 7) = / [(V 5 (r) - 7) - (Vs(r) - /x)] dr 
Jo 

^ — — M7) - KM] 

for M = 2||g||i, and 7 and /1 as in the hypothesis. 

In a similar fashion, let r > 1. Then (z — w) r < 2 r ~ 1 ([[,z] r — Jw] r ) for all 
—M < w < z < M. Indeed, if < w < z, a straightforward argument shows that 

> z T ~ x > (z-wy- 1 ; 

.6 — OJ ^ — UJ 

if w < z < 0, 



and if w < < z, 



> M > (z- wY 

z — w 



■ HI - MT • z r ~\w\ r 1 

mm — ; r — = mm ■ — 



T>6 (z~w) r T>6(z+\w\) r 2 r " 1 
achieved when w = — z. Thus 

(/x - 7 )' r = /' [(Vg(r) - 7) - (Vff(r) - M )f dr 







< T- 1 [^(7) - ■ 

□ 



Theorem 5. Let u and u be solutions of ([T]) with corresponding sources g and g. 
Let m = max{||g||i, ||<?||i}. Then 

\\u - u\\ < 2'- r (\\g - g\\ + (4 "^HIs -g\\l) < r < 1 



||u - fi|| < r2 r - 1 m r - 1 (\\g - g\\ + 2 2 -^ r m 1 - 1 / r r 1 / r \\g - g\\{ /r ^ 



r > 1. 



Proof. Let < r < 1 and s = 2/r > 2. By virtue of ([10 1 

||«-fi|| < ||^([^-7o(g)] r -[^"7o(ff)F)|| s < \\lVg-jo(gW-iv~g-io(~gW\ 

Note that \{zj r - {w} r \ < 2 1 ~ r \z - w\ r for all < \w\ < \z\. Thus 

Hl^ff " lo(gW ~ IV g - 7o(ff)]1 s 1/r < 2^ (\\V(g - g)\\ + \jo(g) - 7o(<?)|) 

<2^ r (||.g-.g|| + |7 (5)- 7o ( ff )|). 
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According to lemma [2J 

ho(~g) - 7o(ff)l < l^(7o(5)) - h g ( l0 (g))\ 



(2m) 1 " 



(2m) 



< 



(4m) 1 -'' „ 

< — — lis- 



[^-7o(5)] r -I^-7o 
(4m) 1 ~ r 



1^(70(9)) -h- g (^(g))\ 



dr 



g\\ r r < 



~g\\ r v 



r r 
This ensures the first statement. 

Let r > 1. In a similar fashion as before, we see that 

||«-fi|| < ||I^ 5 -7o(.9)] r -[^-7o(5)]1 

< r2 r ' 1 m r - 1 (\\V(g - g)\\ + \ l0 (g) - lo (g)\) 

<r2 r - 1 m r - 1 (\\g-g\\ + \ lQ (g)-lo(~g)\)- 
Lemma [4] and similar arguments as for the previous case, yield 

bto® - 7o(ff)l r < r- 1 \h- g { l0 {g)) - h g ( l0 (g))\ 

<2 2r - 2 m r - 1 r\\ g - g \\ 1 . 

This completes the proof. 



□ 



The right hand side bound in the above theorem approaches 2 as r — > indepen- 
dently of the value of \\g — g\\. In fact \\u — u\\ < ||Mli f° r K — Wg~ 7o(ff)] r ~ W(3~ 
7o(g)] r . Since l r (r) -> for almost all r £ [0,1] and |Z r (r)| < 2 + \\Vg\\ 00 + \\Vg\\ 0C , 
the Dominated Convergence Theorem yields \\u — {t|| — > as r — > 0. This is a well- 
known property of the oo-Laplacian, see |22j . As mentioned in the introduction, 
the approach considered in in the context of finite element approximation of 
the solutions of ([!]), may provide an insight on whether the constants found in the 
above theorem are optimal. 



(a) 




10 

q - basis 



2 10 



(b) 




Figure 5. Solving the p-Laplacian problem with p = 5. (a) The 
most accurate basis for a solution u(x) — sin(7ra;) is the standard 
2-sine basis, (b) However for a solution u(x) = fi(x) with q = 5 
the 5-sine basis is the most accurate. 



We now describe how the p-Poisson problem is discretized. The strong formula- 
tion ([T]) leads to the following weak formulation using integration by parts and the 
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boundary conditions: 

(11) f |w'| p "Vv'd:r = f gvdx 

Jo Jo 

for any absolutely continuous test function v. We expand both u and v using a basis 
{4> n }n£Th where <p n is either /„ or /*. After truncation this leads to the following 
nonlinear system of equations for unknown coefficients Cj = (u, <j)*)\ 

p-2 

M'j^ = (9>^) j = l:N. 

The case p = 2 evidently reduces to the standard linear system to solve for the 
Poisson equation. 




(a) (b) 




10 1 10 Z 10 1 

q - basis q - basis 



Figure 6. We now examine a piecewise constant / for different p 
(N = 40). In (a) p = 1.8 and g Q pt = 2 (b) p = 3 and (fopt ~ 2 with 
another minimum at q = 2.95, in (c) p = 5 and q opt ~ 3.9 and in 
(d) p — 10 and q op t ~ 5.75. 



Numerically the right hand side of (12) is approximated via quadrature rules. 
The nonlinear system may then be solved for example by Newton's method. Since 
the Jacobian is a full matrix in this banded approximation appears to 

provide sufficient accuracy. However, the results presented below were found using 
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(a) (b) 




10 1 10 1 
q - basis q - basis 



Figure 7. Comparison of N = 20, 40, 60, 80 in (a) for p = 5 and 
g = 1 (b) for p = 10 and g = 1 (also including N — 100). In (c) 
we show for g = g^ and p = 5, and in (d) for g = g^ and p = 10. 
We note that for N sufficiently large q opt = 2, and a more accurate 
solution may be found for smaller N at q op t ^ 2. 



the trust-region dogleg method with a full Jacobian matrix as implemented in 
Matlab. 

In figure [5] we have use this scheme to solve the p-Laplacian problem with p — 5 
and examine the L 2 error taking N = 40. In figure [El (a) we have fixed a solution 
it(:c) = sin(7rx) where we clearly expect and observe that q = 2 is the optimal basis. 
In figure [51(b) we have fixed u(x) = fi(x) for p — 5 and so we observe the q — 5 as 
the optimal basis. 

In figure[6]we take g — gb from ([£]) which turns out to be a typical form of forcing 
for sandpile problems. We solve the p-Laplacian problem for (a) p — 1.8, (b) p = 3, 
(c) p = 5 and p = 10 with N = 40 and examine the L 2 error in the solution as we 
vary the q basis. To estimate this error we take as exact the solution with 2-sines 
and 2N modes. We observe that the optimal basis for representing the solutions is 
no longer the standard q — 2 for p — 3, 5, 10. For problems with moderate p (for 
example p = 3) we see two distinct minima. For larger p problems however, the 
q = 2 basis becomes less competitive. 

In table [2] we give estimates of q op t for g — gh and g — 1 for N — 20, 40, 60 and 
N = 80 modes. The changes in <7 op t with N are explained by the interchange of the 
two minima in figure [6l For N — 100 this interchange occurs for p = 5 and p = 10 
as illustrated in figure|7l(b). 
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g 


p= 1.8 


p = 'S 


p = 5 


p = 10 


<?b (N = 20) 


2.05 


2.95 


3.9 


8.0 


g h (N = 40) 


2.05 


2.0 (2.95) 


3.9 


5.75 


ffb (N = 60) 


2.2 


3.0 


4.0 


5.95 


<?b (N = 80) 


2.05 


2.0 


3.9 


5.60 


g = l(N = 20) 


2.0 


2.3 


3.1 


6.43 


.g = 1 (AT = 40) 


2.0 


2.0 


3.0 


4.85 


g = l(N = 60) 


2.0 


2.0 


3.0 


4.7 


g = l(N = 80) 


2.0 


2.0 


2.0 (2.95) 


4.7 



Table 2. Optimal q opt for four different p-Laplacian problems 
with either a g — gy, from (|9j) or g = 1. We give estimates for 
different values of N. The change in q opt as iV increases occurs as 
the two minima interchange, see figures [6] and [7j 



5. The time dependent p-Laplacian 

In this final section we consider the evolution equation ^ both in the deter- 
ministic [y = 0) and the stochastically forced regime (i/ ^ 0). Slow- fast diffusion 
is often taken in the large p limit, as a model for sandpile growth. For further 
details including arguments about the validity of the modeling see for example 
[TBI [H [TBI HH HJ [5] • Our purpose here is to consider the time dependent problem 
for different values of p and examine the choice of basis in the formation of the 
sandpile. 

We discretize the weak form of the equation and truncate to solve the time- 



dependent version of ( 12 1 given by 



(13) 



dt 



N 

E 

fc=i 



Cfc 



N 

E 

1=1 



C((j 



p-2 



(b'k^dx - (g, fy) + v(^A 3 ) j = 1 : N. 



This expression is then discretized in time to get a nonlinear system of equations 
to solve for c™ at each step. Below we consider an Euler's method which reduces to 



„n+l 



At 



N 

E' 

fc=i 



N 

E 

i=i 



,n+l. 



P-2 



M 

E 



bj dx 



1/2 



(x)A/3„ 



where A/3 m are independent identically distributed random variables with mean 
zero and variance At (recall ([3])), and we have assumed the eigenfunctions of Q are 
now given by e n . For the case of stochastic forcing we take M » N. 

In figure [8] we consider v = with g = g^ and we plot the time evolution for 
(a) p - 1.8 (<fo p t = 2), (b) p = 3 (cfc p t = 2.95), (c) p = 5 (<z opt = 3.9) and 
(d) p = 10 (<7opt = 5.75). In each case the evolution found numerically quickly 
converges towards the steady state. 

Figure [8] was obtained by choosing in each case the predicted q = q op t basis 
for the steady state. Intuitively this choice of basis should be near optimal for 
large t. On the other hand however, there is no reason to presume that it is so 
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(a) (b) 




Figure 8. Solution of ^ with g = g^ for (a) p = 1.8 using q = 2, 
(b) p = 3 using q = 2.95, (c) p = 5 using g = 3.9 and (d) p = 10 
using g = 5.75. As time evolves the solutions approach the steady 
state which, for increasing p, turns out to be close to a hat function. 



(a) (b) 




Figure 9. Residual in time for the numerical solution of We 
choose v = 0, N = 40 and take as a true solution one computed 
with N — 100 and q = 2. We fix p = 10 and compare basis with 
q = 1.8, q = 2, q = q opt and q = 10. Here (a) corresponds to 
9 = 9b (q op t = 5.75) and (b) to g = 1 (q opt = 4.85). We see 
that asymptotically in time g opt has the smallest error. At a few 
points in time during the transient state, the q = 1.8 basis is more 
accurate than the q = 2. 

for small values of t. In figure [9] we examine the spatial error at each step for 
four different bases with N — 40. For comparison we have chosen p = 10 and 
q G {1-8, 2, q opt , 10}. As the solution approaches the steady state, in each case the 
basis q pt has the smallest error as is expected (in fact a gain of almost an order of 
magnitude in accuracy compared to q = 2 is observed). It is remarkable however 
that, in the transient regime, these bases appear to be far less accurate than other 
choices of q, such as q = 2 and q — 1.8. 

Let us now consider the stochastically forced case, v > 0. Figure[l0]shows plotted 
solution with time dependent noises which are: H 1 in space and white in time for 
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(d) 






0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



(e) 



(f) 



2 0.4 0.6 0.8 1 





Figure 10. Sample realization with p = 10 of the time evolution 
problem with a time-dependent stochastic forcing v = 0.2. The 
time evolution is shown with equal steps of At = 0.01. In (a)-(c) 
we show the effect of an H 1 noise in space, white in time. In (d)-(f) 
we show the effect of a white noise in both space and time. 



(a)-(c), and white both in space and time for (d)-(f). The forcing corresponds 
to v = 0.2. We expect that on average the noisy solution is simply that of the 
deterministic system (which has a unique stable solution). We see in (a)-(b) and 
(d)-(e) the time evolution of the solution for one particular realization of the noise. 
This should be compared with the deterministic case in figure [8j (d). 

In figure 10 -(c) and figure 10 -(f) we plot the evolution of the error for fixed time- 
step with q £ {1, 2, q op t, 10}. Similar to the deterministic case, in the transient 
regime the error for <7 opt and q = 10 is far higher than that for q — 2 or q = 1.8. 
Where the noise effects are large (such as for the spatially correlated noise in H 1 



in figure 10 (c)), the optimal basis is no longer clear. However where the effect is 
small (such as the white noise in figure[lO[(f)), we clearly see in the time dependent 
evolution about the steady state that the q — 4.85 outperforms the other bases (and 
in particular q = 2). For sandpile-type problems the introduction of a spatially 
white noise is a natural choice. It is interesting to note the small effect on the 
dynamics in this realization. 

At present it is unclear whether the use of a g-sine basis with q =/= 2 provides any 
real computational advantage over the natural choice q = 2 for the solution of 
There is clearly a computational overhead in obtaining the former for q ^ 2 that 
impacts significantly on efficiency. However we have observed that the nonlinear 
problem (12 1 is solved faster in the optimal basis. This is certainly worth further 
investigation. For example for TV = 40 we observe approximately a 20% speed up 
in the nonlinear solve over the standard 2-sine basis. Thus, if solving many fixed 
p-problems, the corresponding optimal g-sine basis may not only be more accurate 
but also more efficient. In our current implementation the g-sinc basis may be 
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precomputed and stored. Lemma [T] give some prelimeinary indication on how to 
solve the problem of apriori determining the optimal basis. Our numerical results 
suggest that for large p we expect an optimal basis with q > 2 for a right hand side 
with a discontinuous derivative. 
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